function [ssT_draw,rho_inf, sx , sp,nuhat] = smooth(para,data)

 global psel 
 
 if psel==2 
 
alpha1 = para(1);
alpha2 = para(2);
gamma1 = para(3);
gamma2 = para(4);
beta1  = para(5);
kappa1 = para(6);
tau1   = para(7);
rho_a = para(8);
rho_f = para(9);
rho_i  = para(10);
sd_a1  = para(11);
sd_f1  = para(12);
sd_i1   = para(13);
sd_a2  = para(14);
sd_f2  = para(15);
sd_i2   = para(16);
piess  = para(17);
gy     = para(18);
lnA_0  = para(19);
yst    = para(20);
p11 = para(21);
p22 = para(22);
p11_f = para(23);
p22_f = para(24);

P_i = [p11 1-p11;
       1-p22 p22];
   
P_f = [p11_f 1-p11_f;
       1-p22_f p22_f];   

P = kron(P_i,P_f);

iss    = gy-log(beta1) + piess;

%----------------------------------------------------------------------
% check the roots of the system using gensys
%----------------------------------------------------------------------
A = [p11 1-p11 (p11/tau1) ((1-p11)/tau1);
     1-p22 p22 ((1-p22)/tau1) (p22/tau1);
     0 0 beta1*p11 beta1*(1-p11);
     0 0 beta1*(1-p22) beta1*p22];

B = [1+(gamma1/tau1) 0 (alpha1/tau1) 0;
     0 1+(gamma2/tau1) 0 (alpha2/tau1);
     -kappa1 0 1 0;
     0 -kappa1 0 1];

%----------------------------------------------------------------------
% compute the solutions for inflation and output
%----------------------------------------------------------------------
A_Aa_mat = [1-p11*rho_a, -(1-p11)*rho_a, -(p11*rho_a)/tau1, -((1-p11)*rho_a)/tau1, 1/tau1, 0;
            -(1-p22)*rho_a, 1-p22*rho_a, -((1-p22)*rho_a)/tau1, -(p22*rho_a)/tau1, 0, 1/tau1;
            -kappa1, 0, 1-beta1*p11*rho_a, -beta1*(1-p11)*rho_a, 0, 0;
            0, -kappa1, -beta1*(1-p22)*rho_a, 1-beta1*p22*rho_a, 0, 0;                      
            -gamma1, 0, -alpha1, 0, 1, 0;
            0, -gamma2, 0, -alpha2, 0, 1];           
solvec_a = inv(A_Aa_mat) * [rho_a/tau1;rho_a/tau1;0;0;0;0];

A_Af_mat = [1-p11*rho_f, -(1-p11)*rho_f, -(p11*rho_f)/tau1, -((1-p11)*rho_f)/tau1, 1/tau1, 0;
            -(1-p22)*rho_f, 1-p22*rho_f, -((1-p22)*rho_f)/tau1, -(p22*rho_f)/tau1, 0, 1/tau1;
            -kappa1, 0, 1-beta1*p11*rho_f, -beta1*(1-p11)*rho_f, 0, 0;
            0, -kappa1, -beta1*(1-p22)*rho_f, 1-beta1*p22*rho_f, 0, 0;                      
            -gamma1, 0, -alpha1, 0, 1, 0;
            0, -gamma2, 0, -alpha2, 0, 1];         
solvec_f = inv(A_Af_mat) * [0;0;kappa1;kappa1;0;0];

A_Ai_mat = [1-p11*rho_i, -(1-p11)*rho_i, -(p11*rho_i)/tau1, -((1-p11)*rho_i)/tau1, 1/tau1, 0;
            -(1-p22)*rho_i, 1-p22*rho_i, -((1-p22)*rho_i)/tau1, -(p22*rho_i)/tau1, 0, 1/tau1;
            -kappa1, 0, 1-beta1*p11*rho_i, -beta1*(1-p11)*rho_i, 0, 0;
            0, -kappa1, -beta1*(1-p22)*rho_i, 1-beta1*p22*rho_i, 0, 0;                      
            -gamma1, 0, -alpha1, 0, 1, 0;
            0, -gamma2, 0, -alpha2, 0, 1];        
solvec_i = inv(A_Ai_mat) * [0;0;0;0;1;1];
 
y1a = solvec_a(1);
 y2a = solvec_a(2);
 pi1a = solvec_a(3);
 pi2a = solvec_a(4);
 i1a = solvec_a(5);
 i2a = solvec_a(6);
 y1f = solvec_f(1);
 y2f = solvec_f(2);
 pi1f = solvec_f(3);
 pi2f = solvec_f(4);
 i1f = solvec_f(5);
 i2f = solvec_f(6);
 y1i = solvec_i(1);
 y2i = solvec_i(2);
 pi1i = solvec_i(3);
 pi2i = solvec_i(4);
 i1i=solvec_i(5);
 i2i=solvec_i(6);

%----------------------------------------------------------------------
% compute coefficients for measurement equations 
%----------------------------------------------------------------------
rho     = diag([rho_a;rho_f;rho_i]);

sd1       = [sd_a1;sd_f1;sd_i1];
sd2       = [sd_a2;sd_f2;sd_i2]; 
sd3 =sd1; sd4 =sd2; 


B1 = zeros(1,3); B2 = zeros(1,3);
B3 = zeros(1,3); B4 = zeros(1,3); 

B1(1,:) =  [i1a, i1f, i1i];
B2(1,:) =  [i1a, i1f, i1i];
B3(1,:) =  [i2a, i2f, i2i];
B4(1,:) =  [i2a, i2f, i2i];

% Augment the state vector 
rho1 = diag([diag(rho);1]);
rho2 = diag([diag(rho);1]);
rho3 = diag([diag(rho);1]);
rho4 = diag([diag(rho);1]);
rho1(4,1) = rho_a; 
rho2(4,1) = rho_a; 
rho3(4,1) = rho_a; 
rho4(4,1) = rho_a; 
TC  = [0;0;0;gy]; 
sd1  = [sd1;sd_a1]; sd2  = [sd2;sd_a2];  
sd3  = [sd3;sd_a1]; sd4  = [sd4;sd_a2]; 
Vs1 = diag([sd_a1^2/(1-rho_a^2);sd_f1^2/(1-rho_f^2);sd_i1^2/(1-rho_i^2);0]);
Vs2 = diag([sd_a2^2/(1-rho_a^2);sd_f2^2/(1-rho_f^2);sd_i2^2/(1-rho_i^2);0]);
Vs3 = diag([sd_a1^2/(1-rho_a^2);sd_f1^2/(1-rho_f^2);sd_i1^2/(1-rho_i^2);0]);
Vs4 = diag([sd_a2^2/(1-rho_a^2);sd_f2^2/(1-rho_f^2);sd_i2^2/(1-rho_i^2);0]);

 % Inflation Persistence 
sumw1      = pi1a^2*Vs1(1,1)+pi1f^2*Vs1(2,2)+pi1i^2*Vs1(3,3); 
sumw2      = pi1a^2*Vs2(1,1)+pi1f^2*Vs2(2,2)+pi1i^2*Vs2(3,3);
sumw3      = pi2a^2*Vs3(1,1)+pi2f^2*Vs3(2,2)+pi2i^2*Vs3(3,3); 
sumw4      = pi2a^2*Vs4(1,1)+pi2f^2*Vs4(2,2)+pi2i^2*Vs4(3,3);

weight11   = (pi1a^2)*Vs1(1,1)/sumw1; 
weight12   = (pi1f^2)*Vs1(2,2)/sumw1; 
weight21   = (pi1a^2)*Vs2(1,1)/sumw2; 
weight22   = (pi1f^2)*Vs2(2,2)/sumw2; 
weight31   = (pi2a^2)*Vs3(1,1)/sumw3; 
weight32   = (pi2f^2)*Vs3(2,2)/sumw3; 
weight41   = (pi2a^2)*Vs4(1,1)/sumw4; 
weight42   = (pi2f^2)*Vs4(2,2)/sumw4; 


rho_inf1   = rho_a*weight11+rho_f*weight12+rho_i*(1-weight11-weight12);  
rho_inf2   = rho_a*weight21+rho_f*weight22+rho_i*(1-weight21-weight22); 
rho_inf3   = rho_a*weight31+rho_f*weight32+rho_i*(1-weight31-weight32);  
rho_inf4   = rho_a*weight41+rho_f*weight42+rho_i*(1-weight41-weight42); 


rho_inf    = [rho_inf1;rho_inf2;rho_inf3;rho_inf4]; 

% system matrices for the macro part 
mac1 = [y1a y1f y1i 1; pi1a pi1f pi1i 0]; 
mac2 = [y1a y1f y1i 1; pi1a pi1f pi1i 0]; 
mac3 = [y2a y2f y2i 1; pi2a pi2f pi2i 0]; 
mac4 = [y2a y2f y2i 1; pi2a pi2f pi2i 0]; 

% Markov Chain transition probability 
trans = P';  

% compute the stationary distribution 
A = [eye(length(trans))-trans;ones(1,length(trans))]; 
probs = inv(A'*A)*(A'*[zeros(length(P),1);1]); 
Vs  = probs(1)*Vs1+probs(2)*Vs2+probs(3)*Vs3+probs(4)*Vs4; 

% auxiliary parameters 
nobs = size(data,1);
ns   = size(rho1,1);   % dimension of state vector
nz   = size(data,2);      % dimension of observed variables 
ntr  = ns*(ns+1)/2;
nuhat = zeros(nobs,nz); 
nhist = 2;              
nxmax = 4^nhist; 

Az      = [ yst ; piess ; iss]; 
Az1 = Az; Az2 = Az; Az3 = Az; Az4 = Az; 
Bz1     = [ mac1 ; B1(1,:) 0];    
Bz2     = [ mac2 ; B2(1,:) 0]; 
Bz3     = [ mac3 ; B3(1,:) 0];    
Bz4     = [ mac4 ; B4(1,:) 0]; 


% The matrix for storing information from the forward simulation 
probsT = zeros(nobs,4); 
ssT_draw = zeros(nobs,4); 
probsT(1,:) = probs';
xij     = zeros(16,ns,nobs-1);
xi      = zeros(4,ns,nobs); 
pij     = zeros(16,ntr,nobs-1);
ppi      = zeros(4,ntr,nobs); 

% condition on the first observation : assume the first regime 
s00      = zeros(ns,1); 
s00(4)   = lnA_0; 

llik_c   = zeros(nobs,1);       % log likelihood contribution
nu        = data(1,:)'-(Az+Bz1*(TC+rho1*s00)); 
Ft        = Bz1*Vs*Bz1'; 
llik_c(1) =  -0.5*nz*log(2*pi)-0.5*log(det(Ft))-0.5*nu'*inv(Ft)*nu;

s00      = s00+(Vs*Bz1')*inv(Ft)*nu; 
ltrvecP1 = (ltrvec(Vs-Vs*Bz1'*inv(Ft)*(Vs*Bz1')'))';
probx    = 1;
s1t = s00'; 
xi(1,:,1) = s1t; xi(2,:,1)=s1t; xi(3,:,1)=s1t; xi(4,:,1)=s1t; 
xij(1,:,1) = (TC+rho1*s1t'); xij(2,:,1) = xij(1,:,1); xij(3,:,1)=xij(1,:,1); xij(4,:,1)=xij(1,:,1); 
xij(5,:,1) = xij(1,:,1); xij(6,:,1)=xij(1,:,1); xij(7,:,1)=xij(1,:,1); xij(8,:,1)=xij(1,:,1);  
xij(9,:,1) = xij(1,:,1); xij(10,:,1)=xij(1,:,1); xij(11,:,1)=xij(1,:,1); xij(12,:,1)=xij(1,:,1); 
xij(13,:,1) = xij(1,:,1); xij(14,:,1)=xij(1,:,1); xij(15,:,1)=xij(1,:,1); xij(16,:,1)=xij(1,:,1); 
ppi(1,:,1) = ltrvecP1; ppi(2,:,1)=ppi(1,:,1); ppi(3,:,1) = ppi(1,:,1); ppi(4,:,1)=ppi(1,:,1); 
pij(1,:,1) =  (ltrvec(rho1*ltrmat(ppi(1,:,1)')*rho1'+(diag(sd1))^2))'; pij(2,:,1) =  (ltrvec(rho1*ltrmat(ppi(1,:,1)')*rho1'+(diag(sd2))^2))';
pij(3,:,1) =  (ltrvec(rho1*ltrmat(ppi(1,:,1)')*rho1'+(diag(sd3))^2))'; pij(4,:,1) =  (ltrvec(rho1*ltrmat(ppi(1,:,1)')*rho1'+(diag(sd4))^2))';
pij(5,:,1) =  (ltrvec(rho1*ltrmat(ppi(2,:,1)')*rho1'+(diag(sd1))^2))'; pij(6,:,1) =  (ltrvec(rho1*ltrmat(ppi(2,:,1)')*rho1'+(diag(sd2))^2))';
pij(7,:,1) =  (ltrvec(rho1*ltrmat(ppi(2,:,1)')*rho1'+(diag(sd3))^2))'; pij(8,:,1) =  (ltrvec(rho1*ltrmat(ppi(2,:,1)')*rho1'+(diag(sd4))^2))';
pij(9,:,1) =  (ltrvec(rho1*ltrmat(ppi(3,:,1)')*rho1'+(diag(sd1))^2))'; pij(10,:,1) =  (ltrvec(rho1*ltrmat(ppi(3,:,1)')*rho1'+(diag(sd2))^2))';
pij(11,:,1) = (ltrvec(rho1*ltrmat(ppi(3,:,1)')*rho1'+(diag(sd3))^2))'; pij(12,:,1) =  (ltrvec(rho1*ltrmat(ppi(3,:,1)')*rho1'+(diag(sd4))^2))';
pij(13,:,1) = (ltrvec(rho1*ltrmat(ppi(4,:,1)')*rho1'+(diag(sd1))^2))'; pij(14,:,1) =  (ltrvec(rho1*ltrmat(ppi(4,:,1)')*rho1'+(diag(sd2))^2))';
pij(15,:,1) = (ltrvec(rho1*ltrmat(ppi(4,:,1)')*rho1'+(diag(sd3))^2))'; pij(16,:,1) =  (ltrvec(rho1*ltrmat(ppi(4,:,1)')*rho1'+(diag(sd4))^2))';

t=2;
% log likelihood by Kalman filter
while t<nobs+1
nx    = size(s1t,1);   
probs = trans*probs; 
probx2 = kron(probs,probx); 

llik_ht  = zeros(4*nx,1);      
s0t      = zeros(4*nx,ns);    
ltrvecP0 = zeros(4*nx,ntr);   
yhat     = zeros(nz,1); 
Ft       = zeros(nz,nz); 

i = 1; 
 while i<=nx; 

     
   % forecasting
   
   alphahat1= TC+rho1*s1t(i,:)';
   alphahat2= TC+rho2*s1t(i,:)';
   alphahat3= TC+rho3*s1t(i,:)';
   alphahat4= TC+rho4*s1t(i,:)';
   P1t      = ltrmat((ltrvecP1(i,:))');
   Phat1    = rho1*P1t*rho1'+(diag(sd1))^2; 
   Phat2    = rho2*P1t*rho2'+(diag(sd2))^2; 
   Phat3    = rho3*P1t*rho3'+(diag(sd3))^2; 
   Phat4    = rho4*P1t*rho4'+(diag(sd4))^2; 
   
   yhat1    = Az1 + Bz1*alphahat1;
   yhat2    = Az2 + Bz2*alphahat2; 
   yhat3    = Az3 + Bz3*alphahat3; 
   yhat4    = Az4 + Bz4*alphahat4; 
   
   
   nu1      = data(t,:)-yhat1'; 
   nu2      = data(t,:)-yhat2'; 
   nu3      = data(t,:)-yhat3'; 
   nu4      = data(t,:)-yhat4'; 
   
   Ft1      = Bz1*Phat1*Bz1';
   Ft2      = Bz2*Phat2*Bz2'; 
   Ft3      = Bz3*Phat3*Bz3'; 
   Ft4      = Bz4*Phat4*Bz4'; 
   
   yhat     = yhat + probx2(i)*yhat1+probx2(nx+i)*yhat2+probx2(2*nx+i)*yhat3+probx2(3*nx+i)*yhat4;  
   Ft       = Ft+probx2(i)*(Ft1+yhat1*yhat1')+probx2(nx+i)*(Ft2+yhat2*yhat2')+probx2(2*nx+i)*(Ft3+yhat3*yhat3')+probx2(3*nx+1)*(Ft4+yhat4*yhat4'); 
   
   % compute likelihood 
   
   llik_ht(i,1)    = -0.5*nz*log(2*pi)-0.5*log(det(Ft1))-0.5*nu1*inv(Ft1)*nu1'; 
   llik_ht(nx+i,1) = -0.5*nz*log(2*pi)-0.5*log(det(Ft2))-0.5*nu2*inv(Ft2)*nu2'; 
   llik_ht(2*nx+i,1) = -0.5*nz*log(2*pi)-0.5*log(det(Ft3))-0.5*nu3*inv(Ft3)*nu3'; 
   llik_ht(3*nx+i,1) = -0.5*nz*log(2*pi)-0.5*log(det(Ft4))-0.5*nu4*inv(Ft4)*nu4'; 
    
   % updating step 
   s0t(i,:)    = (alphahat1+Phat1*Bz1'*inv(Ft1)*nu1')';
   s0t(nx+i,:) = (alphahat2+Phat2*Bz2'*inv(Ft2)*nu2')'; 
   s0t(2*nx+i,:) = (alphahat3+Phat3*Bz3'*inv(Ft3)*nu3')'; 
   s0t(3*nx+i,:) = (alphahat4+Phat4*Bz4'*inv(Ft4)*nu4')'; 
   
   ltrvecP0(i,:)=(ltrvec(Phat1-Phat1*Bz1'*inv(Ft1)*Bz1*Phat1'))'; 
   ltrvecP0(nx+i,:) = (ltrvec(Phat2-Phat2*Bz2'*inv(Ft2)*Bz2*Phat2'))'; 
   ltrvecP0(2*nx+i,:) = (ltrvec(Phat3-Phat3*Bz3'*inv(Ft3)*Bz3*Phat3'))'; 
   ltrvecP0(3*nx+i,:) = (ltrvec(Phat4-Phat4*Bz4'*inv(Ft4)*Bz4*Phat4'))'; 
    
   i=i+1;
   
 end 
   
   Ft = Ft -yhat*yhat'; 
   nuhat(t,:) = data(t,:)-yhat';
   
   % llik_c computation 
   
   llik_htmax = max(llik_ht); 
   llik_c(t) = llik_c(t) + llik_htmax + log(probx2'*exp(llik_ht-llik_htmax)); 
   
   % updating mixture probabilities and the state probabilities 
   
   probx2 = (probx2.*exp(llik_ht-llik_htmax))/(probx2'*exp(llik_ht-llik_htmax));
   probs  = [sum(probx2(1:nx));sum(probx2(nx+1:2*nx));sum(probx2(2*nx+1:3*nx));sum(probx2(3*nx+1:4*nx))]; 
   
   % store for forward simulation 
   
   xi(1,:,t) = (probx2(1:nx)'*s0t(1:nx,:))/probs(1); xi(2,:,t) = (probx2(nx+1:2*nx)'*s0t(nx+1:2*nx,:))/probs(2);  
   xi(3,:,t) = (probx2(2*nx+1:3*nx)'*s0t(2*nx+1:3*nx,:))/probs(3); xi(4,:,t) = (probx2(3*nx+1:4*nx)'*s0t(3*nx+1:4*nx,:))/probs(4);  
   xij(1,:,t) = (TC+rho1*xi(1,:,t)'); xij(2,:,t) = (TC+rho1*xi(1,:,t)')'; xij(3,:,t) = (TC+rho1*xi(1,:,t)'); xij(4,:,t) = (TC+rho1*xi(1,:,t)')'; 
   xij(5,:,t) = (TC+rho2*xi(2,:,t)'); xij(6,:,t) = (TC+rho2*xi(2,:,t)')'; xij(7,:,t) = (TC+rho2*xi(2,:,t)'); xij(8,:,t) = (TC+rho1*xi(2,:,t)')'; 
   xij(9,:,t) = (TC+rho3*xi(3,:,t)'); xij(10,:,t) = (TC+rho3*xi(3,:,t)')'; xij(11,:,t) = (TC+rho3*xi(3,:,t)'); xij(12,:,t) = (TC+rho1*xi(3,:,t)')'; 
   xij(13,:,t) = (TC+rho4*xi(4,:,t)'); xij(14,:,t) = (TC+rho4*xi(4,:,t)')'; xij(15,:,t) = (TC+rho4*xi(4,:,t)'); xij(16,:,t) = (TC+rho1*xi(4,:,t)')'; 
   
   paux1  = (ltrvec((s0t(1,:)-xi(1,:,t))'*(s0t(1,:)-xi(1,:,t))))';
   paux2  = (ltrvec((s0t(nx+1,:)-xi(2,:,t))'*(s0t(nx+1,:)-xi(2,:,t))))';
   paux3  = (ltrvec((s0t(2*nx+1,:)-xi(3,:,t))'*(s0t(2*nx+1,:)-xi(3,:,t))))';
   paux4  = (ltrvec((s0t(3*nx+1,:)-xi(4,:,t))'*(s0t(3*nx+1,:)-xi(4,:,t))))';
   
   if nx>1;
      for jj=2:nx 
       paux1 = [paux1;(ltrvec((s0t(jj,:)-xi(1,:,t))'*(s0t(jj,:)-xi(1,:,t))))'];
       paux2 = [paux2;(ltrvec((s0t(nx+jj,:)-xi(2,:,t))'*(s0t(nx+jj,:)-xi(2,:,t))))'];
       paux3 = [paux3;(ltrvec((s0t(2*nx+jj,:)-xi(1,:,t))'*(s0t(2*nx+jj,:)-xi(3,:,t))))'];
       paux4 = [paux4;(ltrvec((s0t(3*nx+jj,:)-xi(2,:,t))'*(s0t(3*nx+jj,:)-xi(4,:,t))))'];
      end 
   end 
   ppi(1,:,t) = probx2(1:nx)'*(ltrvecP0(1:nx,:)+paux1)/probs(1);  ppi(2,:,t) = probx2(nx+1:2*nx)'*(ltrvecP0(nx+1:2*nx,:)+paux2)/probs(2); 
   ppi(3,:,t) = probx2(2*nx+1:3*nx)'*(ltrvecP0(2*nx+1:3*nx,:)+paux3)/probs(3);  ppi(4,:,t) = probx2(3*nx+1:4*nx)'*(ltrvecP0(3*nx+1:4*nx,:)+paux4)/probs(4); 
 
   
   pij(1,:,t) = (ltrvec(rho1*ltrmat(ppi(1,:,t)')*rho1'+(diag(sd1))^2))'; pij(2,:,t) = (ltrvec(rho1*ltrmat(ppi(1,:,t)')*rho1'+(diag(sd2))^2))'; 
   pij(3,:,t) = (ltrvec(rho1*ltrmat(ppi(1,:,t)')*rho1'+(diag(sd3))^2))'; pij(4,:,t) = (ltrvec(rho1*ltrmat(ppi(1,:,t)')*rho1'+(diag(sd4))^2))';   
   pij(5,:,t) = (ltrvec(rho1*ltrmat(ppi(2,:,t)')*rho1'+(diag(sd1))^2))'; pij(6,:,t) = (ltrvec(rho1*ltrmat(ppi(2,:,t)')*rho1'+(diag(sd2))^2))'; 
   pij(7,:,t) = (ltrvec(rho1*ltrmat(ppi(2,:,t)')*rho1'+(diag(sd3))^2))'; pij(8,:,t) = (ltrvec(rho1*ltrmat(ppi(2,:,t)')*rho1'+(diag(sd4))^2))';   
   pij(9,:,t) = (ltrvec(rho1*ltrmat(ppi(3,:,t)')*rho1'+(diag(sd1))^2))'; pij(10,:,t) = (ltrvec(rho1*ltrmat(ppi(3,:,t)')*rho1'+(diag(sd2))^2))'; 
   pij(11,:,t) = (ltrvec(rho1*ltrmat(ppi(3,:,t)')*rho1'+(diag(sd3))^2))'; pij(12,:,t) = (ltrvec(rho1*ltrmat(ppi(3,:,t)')*rho1'+(diag(sd4))^2))';   
   pij(13,:,t) = (ltrvec(rho1*ltrmat(ppi(4,:,t)')*rho1'+(diag(sd1))^2))'; pij(14,:,t) = (ltrvec(rho1*ltrmat(ppi(4,:,t)')*rho1'+(diag(sd2))^2))'; 
   pij(15,:,t) = (ltrvec(rho1*ltrmat(ppi(4,:,t)')*rho1'+(diag(sd3))^2))'; pij(16,:,t) = (ltrvec(rho1*ltrmat(ppi(4,:,t)')*rho1'+(diag(sd4))^2))';   
   
   % collapse the state vector : Kim and Nelson (1999)'s approximation 
   
   zeroprob = (probx2<=1e-5); 
   
   if sum(zeroprob) >= (4*nx-nxmax)
   % simply delete the zero probability states 
   s1t = s0t(zeroprob==0,:);
   ltrvecP1 = ltrvecP0(zeroprob==0,:); 
   probx = probx2(zeroprob==0); 
   
   else

       % collapse some of the states
   selmat       = kron(eye(nx),[1 1 1 1]); 
   sumzeroprob  = selmat*zeroprob;
   del          = sumzeroprob==4;  
   nx4          = sum(del); 
   
   s1t = zeros(nx-nx4,ns);
   ltrvecP1 = zeros(nx-nx4,ntr); 
   probx = zeros(nx-nx4,1); 
   j=1; 
   
   i=1;
   while i<=nx; 
   
       i2 = (i-1)*4+1; 
       
       if sumzeroprob(i)<4 ; 
          % eliminate histories if all the subhistories have insignificant prob.    
          % collapse density 
            probx(j) = probx2(i2)+probx2(i2+1)+probx2(i2+2)+probx2(i2+3); 
            s1t(j,:) = (probx2(i2)/probx(j))*s0t(i2,:)+(probx2(i2+1)/probx(j))*s0t(i2+1,:)+(probx2(i2+2)/probx(j))*s0t(i2+2,:)+(probx2(i2+3)/probx(j))*s0t(i2+3,:);
            ltrvecP1(j,:)=ltrvec((probx2(i2))/(probx(j))*(ltrmat(ltrvecP0(i2,:)')+s0t(i2,:)'*s0t(i2,:)) +...
                          (probx2(i2+1))/(probx(j))*(ltrmat(ltrvecP0(i2+1,:)')+s0t(i2+1,:)'*s0t(i2+1,:))+...
                          (probx2(i2+2))/(probx(j))*(ltrmat(ltrvecP0(i2+2,:)')+s0t(i2+2,:)'*s0t(i2+2,:))+... 
                          (probx2(i2+3))/(probx(j))*(ltrmat(ltrvecP0(i2+3,:)')+s0t(i2+3,:)'*s0t(i2+3,:))+...
                           -s1t(j,:)'*s1t(j,:))';
            j=j+1;
        
       end
       i=i+1;
   end
   end
   probx = probx./sum(probx);
   nnx = size(s1t,1);
      if nnx>nxmax;
       disp('error : nx is greater than nxmax')
       t=nobs+1; 
      end
     % store information for future use 
   probsT(t,:) = probs'; 
t=t+1;
end

    
% Backward Simulation 
t = nobs;
sx.xij = zeros(16,ns,nobs-1);
sx.xi  = zeros(4,ns,nobs);
sx.x   = zeros(nobs,ns); 
sp.pij = zeros(16,ntr,nobs-1);
sp.pi  = zeros(4,ntr,nobs);

sx.xi(:,:,t) = xi(:,:,t); sx.xij(:,:,t-1)=xij(:,:,t-1); sp.pi=ppi(:,:,t); sp.pij(:,:,t-1)=pij(:,:,t-1);
ssT_draw(t,:) = probsT(t,:);
sx.x(t,:) = ssT_draw(t,:)*[sx.xi(1,:,t);sx.xi(2,:,t);sx.xi(3,:,t);sx.xi(4,:,t)]; 

while t>1 
 t = t-1; 
 probs1 = ssT_draw(t+1,:)*trans(:,1)*probsT(t,1)/(ssT_draw(t+1,:)*trans*probsT(t,:)');
 probs2 = ssT_draw(t+1,:)*trans(:,2)*probsT(t,2)/(ssT_draw(t+1,:)*trans*probsT(t,:)');
 probs3 = ssT_draw(t+1,:)*trans(:,3)*probsT(t,3)/(ssT_draw(t+1,:)*trans*probsT(t,:)');
 ssT_draw(t,:) = [probs1, probs2, probs3, 1-probs1-probs2-probs3]; 
 
 for i=1:4 
     for j=1:4 
     paux  = ltrmat(ppi(i,:,t)')*rho1'*inv(ltrmat(pij(4*(i-1)+j,:,t)')); 
     sx.xij(4*(i-1)+j,:,t)=xi(i,:,t)+(paux*(xi(j,:,t+1)-xij(4*(i-1)+j,:,t))')';
     end  
  probaux = [ssT_draw(t+1,1)*probsT(t,i)*trans(1,i), ssT_draw(t+1,2)*probsT(t,i)*trans(2,i) ssT_draw(t+1,3)*probsT(t,i)*trans(3,i) ssT_draw(t+1,4)*probsT(t,i)*trans(4,i)];
  probaux = probaux./sum(probaux);
  sx.xi(i,:,t) = probaux*sx.xij(4*i-3:4*i,:,t); 
  end
  sx.x(t,:) = ssT_draw(t,:)*[sx.xi(1,:,t);sx.xi(2,:,t);sx.xi(3,:,t);sx.xi(4,:,t)]; 
 
end   

 elseif psel ==5 
     
alpha1 = para(1);
alpha2 = para(2);
gamma1 = para(3);
gamma2 = para(4);
beta1  = para(5);
kappa1 = para(6);
tau1   = para(7);
rho_a = para(8);
rho_f = para(9);
rho_i  = para(10);
sd_a1  = para(11);
sd_f1  = para(12);
sd_i1   = para(13);
sd_a2  = para(14);
sd_f2  = para(15);
sd_i2   = para(16);
gy     = para(17);
lnA_0  = para(18);
yst    = para(19);
p11 = para(20);
p22 = para(21);
p11_f = para(22);
p22_f = para(23);
lnP_0  = para(24);
sd_p1  = para(25); 

P_i = [p11 1-p11;
       1-p22 p22];
   
P_f = [p11_f 1-p11_f;
       1-p22_f p22_f];   

P = kron(P_i,P_f);

iss    = gy-log(beta1);

%----------------------------------------------------------------------
% check the roots of the system using gensys
%----------------------------------------------------------------------
A = [p11 1-p11 (p11/tau1) ((1-p11)/tau1);
     1-p22 p22 ((1-p22)/tau1) (p22/tau1);
     0 0 beta1*p11 beta1*(1-p11);
     0 0 beta1*(1-p22) beta1*p22];

B = [1+(gamma1/tau1) 0 (alpha1/tau1) 0;
     0 1+(gamma2/tau1) 0 (alpha2/tau1);
     -kappa1 0 1 0;
     0 -kappa1 0 1];

%----------------------------------------------------------------------
% compute the solutions for inflation and output
%----------------------------------------------------------------------
A_Aa_mat = [1-p11*rho_a, -(1-p11)*rho_a, -(p11*rho_a)/tau1, -((1-p11)*rho_a)/tau1, 1/tau1, 0;
            -(1-p22)*rho_a, 1-p22*rho_a, -((1-p22)*rho_a)/tau1, -(p22*rho_a)/tau1, 0, 1/tau1;
            -kappa1, 0, 1-beta1*p11*rho_a, -beta1*(1-p11)*rho_a, 0, 0;
            0, -kappa1, -beta1*(1-p22)*rho_a, 1-beta1*p22*rho_a, 0, 0;                      
            -gamma1, 0, -alpha1, 0, 1, 0;
            0, -gamma2, 0, -alpha2, 0, 1];           
solvec_a = inv(A_Aa_mat) * [rho_a/tau1;rho_a/tau1;0;0;0;0];

A_Af_mat = [1-p11*rho_f, -(1-p11)*rho_f, -(p11*rho_f)/tau1, -((1-p11)*rho_f)/tau1, 1/tau1, 0;
            -(1-p22)*rho_f, 1-p22*rho_f, -((1-p22)*rho_f)/tau1, -(p22*rho_f)/tau1, 0, 1/tau1;
            -kappa1, 0, 1-beta1*p11*rho_f, -beta1*(1-p11)*rho_f, 0, 0;
            0, -kappa1, -beta1*(1-p22)*rho_f, 1-beta1*p22*rho_f, 0, 0;                      
            -gamma1, 0, -alpha1, 0, 1, 0;
            0, -gamma2, 0, -alpha2, 0, 1];         
solvec_f = inv(A_Af_mat) * [0;0;kappa1;kappa1;0;0];

A_Ai_mat = [1-p11*rho_i, -(1-p11)*rho_i, -(p11*rho_i)/tau1, -((1-p11)*rho_i)/tau1, 1/tau1, 0;
            -(1-p22)*rho_i, 1-p22*rho_i, -((1-p22)*rho_i)/tau1, -(p22*rho_i)/tau1, 0, 1/tau1;
            -kappa1, 0, 1-beta1*p11*rho_i, -beta1*(1-p11)*rho_i, 0, 0;
            0, -kappa1, -beta1*(1-p22)*rho_i, 1-beta1*p22*rho_i, 0, 0;                      
            -gamma1, 0, -alpha1, 0, 1, 0;
            0, -gamma2, 0, -alpha2, 0, 1];        
solvec_i = inv(A_Ai_mat) * [0;0;0;0;1;1];
 
y1a = solvec_a(1);
 y2a = solvec_a(2);
 pi1a = solvec_a(3);
 pi2a = solvec_a(4);
 i1a = solvec_a(5);
 i2a = solvec_a(6);
 y1f = solvec_f(1);
 y2f = solvec_f(2);
 pi1f = solvec_f(3);
 pi2f = solvec_f(4);
 i1f = solvec_f(5);
 i2f = solvec_f(6);
 y1i = solvec_i(1);
 y2i = solvec_i(2);
 pi1i = solvec_i(3);
 pi2i = solvec_i(4);
 i1i=solvec_i(5);
 i2i=solvec_i(6);

%----------------------------------------------------------------------
% compute coefficients for measurement equations 
%----------------------------------------------------------------------
rho     = diag([rho_a;rho_f;rho_i]);

sd1       = [sd_a1;sd_f1;sd_i1];
sd2       = [sd_a2;sd_f2;sd_i2]; 
sd3 =sd1; sd4 =sd2; 


B1 = zeros(1,3); B2 = zeros(1,3);
B3 = zeros(1,3); B4 = zeros(1,3); 

B1(1,:) =  [i1a, i1f, i1i];
B2(1,:) =  [i1a, i1f, i1i];
B3(1,:) =  [i2a, i2f, i2i];
B4(1,:) =  [i2a, i2f, i2i];

% Augment the state vector 
rho1 = diag([diag(rho);1;1]);
rho2 = diag([diag(rho);1;1]);
rho3 = diag([diag(rho);1;1]);
rho4 = diag([diag(rho);1;1]);
rho1(4,1) = rho_a; 
rho2(4,1) = rho_a; 
rho3(4,1) = rho_a; 
rho4(4,1) = rho_a; 
TC  = [0;0;0;gy;0]; 
sd1  = [sd1;sd_a1;sd_p1]; sd2  = [sd2;sd_a2;sd_p1];  
sd3  = [sd3;sd_a1;sd_p1]; sd4  = [sd4;sd_a2;sd_p1]; 
Vs1 = diag([sd_a1^2/(1-rho_a^2);sd_f1^2/(1-rho_f^2);sd_i1^2/(1-rho_i^2);0;0]);
Vs2 = diag([sd_a2^2/(1-rho_a^2);sd_f2^2/(1-rho_f^2);sd_i2^2/(1-rho_i^2);0;0]);
Vs3 = diag([sd_a1^2/(1-rho_a^2);sd_f1^2/(1-rho_f^2);sd_i1^2/(1-rho_i^2);0;0]);
Vs4 = diag([sd_a2^2/(1-rho_a^2);sd_f2^2/(1-rho_f^2);sd_i2^2/(1-rho_i^2);0;0]);

 % Inflation Persistence 
sumw1      = pi1a^2*Vs1(1,1)+pi1f^2*Vs1(2,2)+pi1i^2*Vs1(3,3); 
sumw2      = pi1a^2*Vs2(1,1)+pi1f^2*Vs2(2,2)+pi1i^2*Vs2(3,3);
sumw3      = pi2a^2*Vs3(1,1)+pi2f^2*Vs3(2,2)+pi2i^2*Vs3(3,3); 
sumw4      = pi2a^2*Vs4(1,1)+pi2f^2*Vs4(2,2)+pi2i^2*Vs4(3,3);

weight11   = (pi1a^2)*Vs1(1,1)/sumw1; 
weight12   = (pi1f^2)*Vs1(2,2)/sumw1; 
weight21   = (pi1a^2)*Vs2(1,1)/sumw2; 
weight22   = (pi1f^2)*Vs2(2,2)/sumw2; 
weight31   = (pi2a^2)*Vs3(1,1)/sumw3; 
weight32   = (pi2f^2)*Vs3(2,2)/sumw3; 
weight41   = (pi2a^2)*Vs4(1,1)/sumw4; 
weight42   = (pi2f^2)*Vs4(2,2)/sumw4; 


rho_inf1   = rho_a*weight11+rho_f*weight12+rho_i*(1-weight11-weight12);  
rho_inf2   = rho_a*weight21+rho_f*weight22+rho_i*(1-weight21-weight22); 
rho_inf3   = rho_a*weight31+rho_f*weight32+rho_i*(1-weight31-weight32);  
rho_inf4   = rho_a*weight41+rho_f*weight42+rho_i*(1-weight41-weight42); 


rho_inf    = [rho_inf1;rho_inf2;rho_inf3;rho_inf4]; 

% system matrices for the macro part 
mac1 = [y1a y1f y1i 1 0; pi1a pi1f pi1i 0 1]; 
mac2 = [y1a y1f y1i 1 0; pi1a pi1f pi1i 0 1]; 
mac3 = [y2a y2f y2i 1 0; pi2a pi2f pi2i 0 1]; 
mac4 = [y2a y2f y2i 1 0; pi2a pi2f pi2i 0 1]; 

% Markov Chain transition probability 
trans = P';  

% compute the stationary distribution 
A = [eye(length(trans))-trans;ones(1,length(trans))]; 
probs = inv(A'*A)*(A'*[zeros(length(P),1);1]); 
Vs  = probs(1)*Vs1+probs(2)*Vs2+probs(3)*Vs3+probs(4)*Vs4; 

% auxiliary parameters 
nobs = size(data,1);
ns   = size(rho1,1);   % dimension of state vector
nz   = size(data,2);      % dimension of observed variables 
ntr  = ns*(ns+1)/2;
nuhat = zeros(nobs,nz); 
nhist = 2;              
nxmax = 4^nhist; 

Az      = [ yst ; 0 ; iss]; 
Az1 = Az; Az2 = Az; Az3 = Az; Az4 = Az; 
Bz1     = [ mac1 ; B1(1,:) 0 1];    
Bz2     = [ mac2 ; B2(1,:) 0 1]; 
Bz3     = [ mac3 ; B3(1,:) 0 1];    
Bz4     = [ mac4 ; B4(1,:) 0 1]; 


% The matrix for storing information from the forward simulation 
probsT = zeros(nobs,4); 
ssT_draw = zeros(nobs,4); 
probsT(1,:) = probs';
xij     = zeros(16,ns,nobs-1);
xi      = zeros(4,ns,nobs); 
pij     = zeros(16,ntr,nobs-1);
ppi      = zeros(4,ntr,nobs); 

% condition on the first observation : assume the first regime 
s00      = zeros(ns,1); 
s00(4)   = lnA_0; 
s00(5)   = lnP_0; 

llik_c   = zeros(nobs,1);       % log likelihood contribution
nu        = data(1,:)'-(Az+Bz1*(TC+rho1*s00)); 
Ft        = Bz1*Vs*Bz1'; 
llik_c(1) =  -0.5*nz*log(2*pi)-0.5*log(det(Ft))-0.5*nu'*inv(Ft)*nu;

s00      = s00+(Vs*Bz1')*inv(Ft)*nu; 
ltrvecP1 = (ltrvec(Vs-Vs*Bz1'*inv(Ft)*(Vs*Bz1')'))';
probx    = 1;
s1t = s00'; 
xi(1,:,1) = s1t; xi(2,:,1)=s1t; xi(3,:,1)=s1t; xi(4,:,1)=s1t; 
xij(1,:,1) = (TC+rho1*s1t'); xij(2,:,1) = xij(1,:,1); xij(3,:,1)=xij(1,:,1); xij(4,:,1)=xij(1,:,1); 
xij(5,:,1) = xij(1,:,1); xij(6,:,1)=xij(1,:,1); xij(7,:,1)=xij(1,:,1); xij(8,:,1)=xij(1,:,1);  
xij(9,:,1) = xij(1,:,1); xij(10,:,1)=xij(1,:,1); xij(11,:,1)=xij(1,:,1); xij(12,:,1)=xij(1,:,1); 
xij(13,:,1) = xij(1,:,1); xij(14,:,1)=xij(1,:,1); xij(15,:,1)=xij(1,:,1); xij(16,:,1)=xij(1,:,1); 
ppi(1,:,1) = ltrvecP1; ppi(2,:,1)=ppi(1,:,1); ppi(3,:,1) = ppi(1,:,1); ppi(4,:,1)=ppi(1,:,1); 
pij(1,:,1) =  (ltrvec(rho1*ltrmat(ppi(1,:,1)')*rho1'+(diag(sd1))^2))'; pij(2,:,1) =  (ltrvec(rho1*ltrmat(ppi(1,:,1)')*rho1'+(diag(sd2))^2))';
pij(3,:,1) =  (ltrvec(rho1*ltrmat(ppi(1,:,1)')*rho1'+(diag(sd3))^2))'; pij(4,:,1) =  (ltrvec(rho1*ltrmat(ppi(1,:,1)')*rho1'+(diag(sd4))^2))';
pij(5,:,1) =  (ltrvec(rho1*ltrmat(ppi(2,:,1)')*rho1'+(diag(sd1))^2))'; pij(6,:,1) =  (ltrvec(rho1*ltrmat(ppi(2,:,1)')*rho1'+(diag(sd2))^2))';
pij(7,:,1) =  (ltrvec(rho1*ltrmat(ppi(2,:,1)')*rho1'+(diag(sd3))^2))'; pij(8,:,1) =  (ltrvec(rho1*ltrmat(ppi(2,:,1)')*rho1'+(diag(sd4))^2))';
pij(9,:,1) =  (ltrvec(rho1*ltrmat(ppi(3,:,1)')*rho1'+(diag(sd1))^2))'; pij(10,:,1) =  (ltrvec(rho1*ltrmat(ppi(3,:,1)')*rho1'+(diag(sd2))^2))';
pij(11,:,1) = (ltrvec(rho1*ltrmat(ppi(3,:,1)')*rho1'+(diag(sd3))^2))'; pij(12,:,1) =  (ltrvec(rho1*ltrmat(ppi(3,:,1)')*rho1'+(diag(sd4))^2))';
pij(13,:,1) = (ltrvec(rho1*ltrmat(ppi(4,:,1)')*rho1'+(diag(sd1))^2))'; pij(14,:,1) =  (ltrvec(rho1*ltrmat(ppi(4,:,1)')*rho1'+(diag(sd2))^2))';
pij(15,:,1) = (ltrvec(rho1*ltrmat(ppi(4,:,1)')*rho1'+(diag(sd3))^2))'; pij(16,:,1) =  (ltrvec(rho1*ltrmat(ppi(4,:,1)')*rho1'+(diag(sd4))^2))';

t=2;
% log likelihood by Kalman filter
while t<nobs+1
nx    = size(s1t,1);   
probs = trans*probs; 
probx2 = kron(probs,probx); 

llik_ht  = zeros(4*nx,1);      
s0t      = zeros(4*nx,ns);    
ltrvecP0 = zeros(4*nx,ntr);   
yhat     = zeros(nz,1); 
Ft       = zeros(nz,nz); 

i = 1; 
 while i<=nx; 

     
   % forecasting
   
   alphahat1= TC+rho1*s1t(i,:)';
   alphahat2= TC+rho2*s1t(i,:)';
   alphahat3= TC+rho3*s1t(i,:)';
   alphahat4= TC+rho4*s1t(i,:)';
   P1t      = ltrmat((ltrvecP1(i,:))');
   Phat1    = rho1*P1t*rho1'+(diag(sd1))^2; 
   Phat2    = rho2*P1t*rho2'+(diag(sd2))^2; 
   Phat3    = rho3*P1t*rho3'+(diag(sd3))^2; 
   Phat4    = rho4*P1t*rho4'+(diag(sd4))^2; 
   
   yhat1    = Az1 + Bz1*alphahat1;
   yhat2    = Az2 + Bz2*alphahat2; 
   yhat3    = Az3 + Bz3*alphahat3; 
   yhat4    = Az4 + Bz4*alphahat4; 
   
   
   nu1      = data(t,:)-yhat1'; 
   nu2      = data(t,:)-yhat2'; 
   nu3      = data(t,:)-yhat3'; 
   nu4      = data(t,:)-yhat4'; 
   
   Ft1      = Bz1*Phat1*Bz1';
   Ft2      = Bz2*Phat2*Bz2'; 
   Ft3      = Bz3*Phat3*Bz3'; 
   Ft4      = Bz4*Phat4*Bz4'; 
   
   yhat     = yhat + probx2(i)*yhat1+probx2(nx+i)*yhat2+probx2(2*nx+i)*yhat3+probx2(3*nx+i)*yhat4;  
   Ft       = Ft+probx2(i)*(Ft1+yhat1*yhat1')+probx2(nx+i)*(Ft2+yhat2*yhat2')+probx2(2*nx+i)*(Ft3+yhat3*yhat3')+probx2(3*nx+1)*(Ft4+yhat4*yhat4'); 
   
   % compute likelihood 
   
   llik_ht(i,1)    = -0.5*nz*log(2*pi)-0.5*log(det(Ft1))-0.5*nu1*inv(Ft1)*nu1'; 
   llik_ht(nx+i,1) = -0.5*nz*log(2*pi)-0.5*log(det(Ft2))-0.5*nu2*inv(Ft2)*nu2'; 
   llik_ht(2*nx+i,1) = -0.5*nz*log(2*pi)-0.5*log(det(Ft3))-0.5*nu3*inv(Ft3)*nu3'; 
   llik_ht(3*nx+i,1) = -0.5*nz*log(2*pi)-0.5*log(det(Ft4))-0.5*nu4*inv(Ft4)*nu4'; 
    
   % updating step 
   s0t(i,:)    = (alphahat1+Phat1*Bz1'*inv(Ft1)*nu1')';
   s0t(nx+i,:) = (alphahat2+Phat2*Bz2'*inv(Ft2)*nu2')'; 
   s0t(2*nx+i,:) = (alphahat3+Phat3*Bz3'*inv(Ft3)*nu3')'; 
   s0t(3*nx+i,:) = (alphahat4+Phat4*Bz4'*inv(Ft4)*nu4')'; 
   
   ltrvecP0(i,:)=(ltrvec(Phat1-Phat1*Bz1'*inv(Ft1)*Bz1*Phat1'))'; 
   ltrvecP0(nx+i,:) = (ltrvec(Phat2-Phat2*Bz2'*inv(Ft2)*Bz2*Phat2'))'; 
   ltrvecP0(2*nx+i,:) = (ltrvec(Phat3-Phat3*Bz3'*inv(Ft3)*Bz3*Phat3'))'; 
   ltrvecP0(3*nx+i,:) = (ltrvec(Phat4-Phat4*Bz4'*inv(Ft4)*Bz4*Phat4'))'; 
    
   i=i+1;
   
 end 
   
   Ft = Ft -yhat*yhat'; 
   nuhat(t,:) = data(t,:)-yhat';
   
   % llik_c computation 
   
   llik_htmax = max(llik_ht); 
   llik_c(t) = llik_c(t) + llik_htmax + log(probx2'*exp(llik_ht-llik_htmax)); 
   
   % updating mixture probabilities and the state probabilities 
   
   probx2 = (probx2.*exp(llik_ht-llik_htmax))/(probx2'*exp(llik_ht-llik_htmax));
   probs  = [sum(probx2(1:nx));sum(probx2(nx+1:2*nx));sum(probx2(2*nx+1:3*nx));sum(probx2(3*nx+1:4*nx))]; 
   
   % store for forward simulation 
   
   xi(1,:,t) = (probx2(1:nx)'*s0t(1:nx,:))/probs(1); xi(2,:,t) = (probx2(nx+1:2*nx)'*s0t(nx+1:2*nx,:))/probs(2);  
   xi(3,:,t) = (probx2(2*nx+1:3*nx)'*s0t(2*nx+1:3*nx,:))/probs(3); xi(4,:,t) = (probx2(3*nx+1:4*nx)'*s0t(3*nx+1:4*nx,:))/probs(4);  
   xij(1,:,t) = (TC+rho1*xi(1,:,t)'); xij(2,:,t) = (TC+rho1*xi(1,:,t)')'; xij(3,:,t) = (TC+rho1*xi(1,:,t)'); xij(4,:,t) = (TC+rho1*xi(1,:,t)')'; 
   xij(5,:,t) = (TC+rho2*xi(2,:,t)'); xij(6,:,t) = (TC+rho2*xi(2,:,t)')'; xij(7,:,t) = (TC+rho2*xi(2,:,t)'); xij(8,:,t) = (TC+rho1*xi(2,:,t)')'; 
   xij(9,:,t) = (TC+rho3*xi(3,:,t)'); xij(10,:,t) = (TC+rho3*xi(3,:,t)')'; xij(11,:,t) = (TC+rho3*xi(3,:,t)'); xij(12,:,t) = (TC+rho1*xi(3,:,t)')'; 
   xij(13,:,t) = (TC+rho4*xi(4,:,t)'); xij(14,:,t) = (TC+rho4*xi(4,:,t)')'; xij(15,:,t) = (TC+rho4*xi(4,:,t)'); xij(16,:,t) = (TC+rho1*xi(4,:,t)')'; 
   
   paux1  = (ltrvec((s0t(1,:)-xi(1,:,t))'*(s0t(1,:)-xi(1,:,t))))';
   paux2  = (ltrvec((s0t(nx+1,:)-xi(2,:,t))'*(s0t(nx+1,:)-xi(2,:,t))))';
   paux3  = (ltrvec((s0t(2*nx+1,:)-xi(3,:,t))'*(s0t(2*nx+1,:)-xi(3,:,t))))';
   paux4  = (ltrvec((s0t(3*nx+1,:)-xi(4,:,t))'*(s0t(3*nx+1,:)-xi(4,:,t))))';
   
   if nx>1;
      for jj=2:nx 
       paux1 = [paux1;(ltrvec((s0t(jj,:)-xi(1,:,t))'*(s0t(jj,:)-xi(1,:,t))))'];
       paux2 = [paux2;(ltrvec((s0t(nx+jj,:)-xi(2,:,t))'*(s0t(nx+jj,:)-xi(2,:,t))))'];
       paux3 = [paux3;(ltrvec((s0t(2*nx+jj,:)-xi(1,:,t))'*(s0t(2*nx+jj,:)-xi(3,:,t))))'];
       paux4 = [paux4;(ltrvec((s0t(3*nx+jj,:)-xi(2,:,t))'*(s0t(3*nx+jj,:)-xi(4,:,t))))'];
      end 
   end 
   ppi(1,:,t) = probx2(1:nx)'*(ltrvecP0(1:nx,:)+paux1)/probs(1);  ppi(2,:,t) = probx2(nx+1:2*nx)'*(ltrvecP0(nx+1:2*nx,:)+paux2)/probs(2); 
   ppi(3,:,t) = probx2(2*nx+1:3*nx)'*(ltrvecP0(2*nx+1:3*nx,:)+paux3)/probs(3);  ppi(4,:,t) = probx2(3*nx+1:4*nx)'*(ltrvecP0(3*nx+1:4*nx,:)+paux4)/probs(4); 
 
   
   pij(1,:,t) = (ltrvec(rho1*ltrmat(ppi(1,:,t)')*rho1'+(diag(sd1))^2))'; pij(2,:,t) = (ltrvec(rho1*ltrmat(ppi(1,:,t)')*rho1'+(diag(sd2))^2))'; 
   pij(3,:,t) = (ltrvec(rho1*ltrmat(ppi(1,:,t)')*rho1'+(diag(sd3))^2))'; pij(4,:,t) = (ltrvec(rho1*ltrmat(ppi(1,:,t)')*rho1'+(diag(sd4))^2))';   
   pij(5,:,t) = (ltrvec(rho1*ltrmat(ppi(2,:,t)')*rho1'+(diag(sd1))^2))'; pij(6,:,t) = (ltrvec(rho1*ltrmat(ppi(2,:,t)')*rho1'+(diag(sd2))^2))'; 
   pij(7,:,t) = (ltrvec(rho1*ltrmat(ppi(2,:,t)')*rho1'+(diag(sd3))^2))'; pij(8,:,t) = (ltrvec(rho1*ltrmat(ppi(2,:,t)')*rho1'+(diag(sd4))^2))';   
   pij(9,:,t) = (ltrvec(rho1*ltrmat(ppi(3,:,t)')*rho1'+(diag(sd1))^2))'; pij(10,:,t) = (ltrvec(rho1*ltrmat(ppi(3,:,t)')*rho1'+(diag(sd2))^2))'; 
   pij(11,:,t) = (ltrvec(rho1*ltrmat(ppi(3,:,t)')*rho1'+(diag(sd3))^2))'; pij(12,:,t) = (ltrvec(rho1*ltrmat(ppi(3,:,t)')*rho1'+(diag(sd4))^2))';   
   pij(13,:,t) = (ltrvec(rho1*ltrmat(ppi(4,:,t)')*rho1'+(diag(sd1))^2))'; pij(14,:,t) = (ltrvec(rho1*ltrmat(ppi(4,:,t)')*rho1'+(diag(sd2))^2))'; 
   pij(15,:,t) = (ltrvec(rho1*ltrmat(ppi(4,:,t)')*rho1'+(diag(sd3))^2))'; pij(16,:,t) = (ltrvec(rho1*ltrmat(ppi(4,:,t)')*rho1'+(diag(sd4))^2))';   
   
   % collapse the state vector : Kim and Nelson (1999)'s approximation 
   
   zeroprob = (probx2<=1e-5); 
   
   if sum(zeroprob) >= (4*nx-nxmax)
   % simply delete the zero probability states 
   s1t = s0t(zeroprob==0,:);
   ltrvecP1 = ltrvecP0(zeroprob==0,:); 
   probx = probx2(zeroprob==0); 
   
   else

       % collapse some of the states
   selmat       = kron(eye(nx),[1 1 1 1]); 
   sumzeroprob  = selmat*zeroprob;
   del          = sumzeroprob==4;  
   nx4          = sum(del); 
   
   s1t = zeros(nx-nx4,ns);
   ltrvecP1 = zeros(nx-nx4,ntr); 
   probx = zeros(nx-nx4,1); 
   j=1; 
   
   i=1;
   while i<=nx; 
   
       i2 = (i-1)*4+1; 
       
       if sumzeroprob(i)<4 ; 
          % eliminate histories if all the subhistories have insignificant prob.    
          % collapse density 
            probx(j) = probx2(i2)+probx2(i2+1)+probx2(i2+2)+probx2(i2+3); 
            s1t(j,:) = (probx2(i2)/probx(j))*s0t(i2,:)+(probx2(i2+1)/probx(j))*s0t(i2+1,:)+(probx2(i2+2)/probx(j))*s0t(i2+2,:)+(probx2(i2+3)/probx(j))*s0t(i2+3,:);
            ltrvecP1(j,:)=ltrvec((probx2(i2))/(probx(j))*(ltrmat(ltrvecP0(i2,:)')+s0t(i2,:)'*s0t(i2,:)) +...
                          (probx2(i2+1))/(probx(j))*(ltrmat(ltrvecP0(i2+1,:)')+s0t(i2+1,:)'*s0t(i2+1,:))+...
                          (probx2(i2+2))/(probx(j))*(ltrmat(ltrvecP0(i2+2,:)')+s0t(i2+2,:)'*s0t(i2+2,:))+... 
                          (probx2(i2+3))/(probx(j))*(ltrmat(ltrvecP0(i2+3,:)')+s0t(i2+3,:)'*s0t(i2+3,:))+...
                           -s1t(j,:)'*s1t(j,:))';
            j=j+1;
        
       end
       i=i+1;
   end
   end
   probx = probx./sum(probx);
   nnx = size(s1t,1);
      if nnx>nxmax;
       disp('error : nx is greater than nxmax')
       t=nobs+1; 
      end
     % store information for future use 
   probsT(t,:) = probs'; 
t=t+1;
end

    
% Backward Simulation 
t = nobs;
sx.xij = zeros(16,ns,nobs-1);
sx.xi  = zeros(4,ns,nobs);
sx.x   = zeros(nobs,ns); 
sp.pij = zeros(16,ntr,nobs-1);
sp.pi  = zeros(4,ntr,nobs);

sx.xi(:,:,t) = xi(:,:,t); sx.xij(:,:,t-1)=xij(:,:,t-1); sp.pi=ppi(:,:,t); sp.pij(:,:,t-1)=pij(:,:,t-1);
ssT_draw(t,:) = probsT(t,:);
sx.x(t,:) = ssT_draw(t,:)*[sx.xi(1,:,t);sx.xi(2,:,t);sx.xi(3,:,t);sx.xi(4,:,t)]; 

while t>1 
 t = t-1; 
 probs1 = ssT_draw(t+1,:)*trans(:,1)*probsT(t,1)/(ssT_draw(t+1,:)*trans*probsT(t,:)');
 probs2 = ssT_draw(t+1,:)*trans(:,2)*probsT(t,2)/(ssT_draw(t+1,:)*trans*probsT(t,:)');
 probs3 = ssT_draw(t+1,:)*trans(:,3)*probsT(t,3)/(ssT_draw(t+1,:)*trans*probsT(t,:)');
 ssT_draw(t,:) = [probs1, probs2, probs3, 1-probs1-probs2-probs3]; 
 
 for i=1:4 
     for j=1:4 
     paux  = ltrmat(ppi(i,:,t)')*rho1'*inv(ltrmat(pij(4*(i-1)+j,:,t)')); 
     sx.xij(4*(i-1)+j,:,t)=xi(i,:,t)+(paux*(xi(j,:,t+1)-xij(4*(i-1)+j,:,t))')';
     end  
  probaux = [ssT_draw(t+1,1)*probsT(t,i)*trans(1,i), ssT_draw(t+1,2)*probsT(t,i)*trans(2,i) ssT_draw(t+1,3)*probsT(t,i)*trans(3,i) ssT_draw(t+1,4)*probsT(t,i)*trans(4,i)];
  probaux = probaux./sum(probaux);
  sx.xi(i,:,t) = probaux*sx.xij(4*i-3:4*i,:,t); 
  end
  sx.x(t,:) = ssT_draw(t,:)*[sx.xi(1,:,t);sx.xi(2,:,t);sx.xi(3,:,t);sx.xi(4,:,t)]; 
 
end   
     
 end 
  